classdef nanSkewness < handle
  % Derived from handle
  % Constructor signature: skewness = nanSkewness(x)
  %
  % Substantive methods:
  %
  %     Pearson2 = Pearson's Modified Second Coefficient of Skewness
  %              = (mean - meadian) / standard deviation
  %         ignoring NaNs
  %    
  %     Fisher = "Fisher's Skewness"
  %            = n * sum(z^3) / [(n-1)(n-2)]
  %         ignoring NaNs
  %
  %     TwoSidedFisherSkewnessTest(a) implements a two-sided version of the
  %         test described in the Reference, following the formula it gives
  %         for Fisher's Skewness; a is the desired confidence level.
  %
  %     Reference: Wuensch, Karl L., 2007. "Skewness, Kurtosis, and the Normal Curve."
  %                http://core.ecu.edu/psyc/wuenschk/docs30/Skew-Kurt.doc
  %
  % Author: David Goldsmith, Wash. State Dept. of Ecology, dgol461@ecy.wa.gov
  % Release date: 4/20/2011    

    properties
        good; n; x; meen; med; s; skew;
    end

    methods
        function skewness = nanSkewness(x)
            skewness.good = ~(isnan(x) | isinf(abs(x)));
            skewness.n    = sum(skewness.good);
            skewness.x    = x(skewness.good);
            skewness.meen = mean(skewness.x);
            skewness.med  = median(skewness.x);
            skewness.s    = std(skewness.x);
        end
        
        function Pearson2(obj)
            if obj.n
                obj.skew = (obj.meen - obj.med) / obj.s;
            else
                obj.skew = NaN;
            end
        end
        
        function Fisher(obj)
            if obj.n < 3
                obj.skew = NaN;
            else
                Sz3 = sum(((obj.x - obj.meen) / obj.s).^3);
                obj.skew = obj.n * Sz3 / ((obj.n-1)*(obj.n-2));
            end
        end
        
        function [isskew, y] = TwoSidedFisherSkewnessTest(obj, a)
            if obj.n <= 150
                y = 0;
                isskew = NaN;
                return
            end
            
            if nargin < 2
                a = 0.95; % Default confidence level
            end
            
          % Compute Fisher Skewness
            obj.Fisher();
            
          % According to Reference, n > 150 implies Fisher Skewness is 
          % approx. normal with standard error of approx. sqrt(6/n);
          % y = Z(a, 2-sided) / that s. e.
            y = sqrt(obj.n/3) * obj.skew / erfinv(a) / 2;
            if abs(y) > 5
                isskew = 1;
            else
                isskew = 0;
            end
        end
    end
end